/*
Este programa implementa o algoritmo de Jacobi-Richardson para a resolução iterativa
de sistemas lineares (http://en.wikipedia.org/wiki/Jacobi_method).
A versão contida neste arquivo é paralelizada em CUDA.
*/

#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>

//Tamanho de cada bloco na grid. É vantajoso, nesse caso, usar o máximo
//disponível, que no caso da placa de vídeo utilizada é 512.
#define BLOCK_SIZE 512

//Tamanho máximo da grid. Este parâmetro não é utilizado no programa, mas
//é útil caso se desejem tratar matrizes de tamanho maior ou igual a este.
#define MAX_GRID_SIZE 65535


//Esta função calcula o produto Ax (com o valor estimado de x na iteração
//atual) e compara com os valores do vetor b. Os erros são guardados no vetor
//err. Cada bloco é responsável pelo produto de uma linha da matriz pelo vetor,
//e uma redução é realizada ao final por uma thread de cada bloco.
__global__ void errors(float *A, float *x, float *b, float *err, int n)
{
   int i, j;
   __shared__ float results[BLOCK_SIZE];
   float thread_result;
   
   j = blockIdx.x;
   thread_result = 0;
   
   for(i = threadIdx.x; i < n; i += BLOCK_SIZE)
      thread_result += A[j * n + i] * x[i];
   results[threadIdx.x] = thread_result;
   __syncthreads();
   
   if(threadIdx.x == 0)
   {
      for(i = 1; i < BLOCK_SIZE; i++)
         thread_result += results[i];
      err[j] = b[j] - thread_result;
   }
}

//Esta função é responsável por atualizar o valor de x ao da iteração atual.
//Utiliza o valor do vetor x_p, da iteração anterior, calculando x = B * x_p + g.
//Cada bloco é responsável pelo cálculo de um elemento do vetor x, e uma redução
//é realizada ao final por uma thread de cada bloco. 
__global__ void update(float *B, float *g, float *x, float *x_p, int n)
{
   int i, j;
   __shared__ float results[BLOCK_SIZE];
   float thread_result;
   
   j = blockIdx.x;
   thread_result = 0;
   
   for(i = threadIdx.x; i < n; i += BLOCK_SIZE)
      thread_result += B[j * n + i] * x_p[i];
   results[threadIdx.x] = thread_result;
   __syncthreads();
   
   if(threadIdx.x == 0)
   {
      for(i = 1; i < BLOCK_SIZE; i++)
         thread_result += results[i];
      x[j] = thread_result + g[j];
   }
}

//Esta função é responsável pelo cálculo das matrizes B e g, a partir de A e b.
//Todos os elementos de cada linha (e do vetor b) são divididos pelo valor do
//elemento na diagonal da respectiva linha, e os elementos da diagonal de B
//são 0. Cada bloco calcula uma linha. Também inicializa o vetor x com zeros.
__global__ void setup(float *A, float *b, float *B, float *g, float *x, int n)
{
   int i, j;
   
   j = blockIdx.x;
   
   for(i = threadIdx.x; i < n; i += BLOCK_SIZE)
      if(i != j)
         B[j * n + i] = -(A[j * n + i] / A[j * n + j]);
      else
         B[j * n + i] = 0;
         
   if(blockIdx.x == 0)
      for(i = threadIdx.x; i < n; i += BLOCK_SIZE)
      {
         g[i] = b[i] / A[i * n + i];
         x[i] = 0;
      }
}

int main(int argc, const char* argv[])
{
   struct timeval ti, tf;
   long int sec, usec;
   //Início da contagem de tempo
   gettimeofday(&ti, NULL);
   
   
   float *hA, *hb, *dA, *db, *herrors;
   float *dB, *dg;
   float *hx, *dx[2], *derrors;
   int n, row, iterations, max_iterations;
   float error, max_error;
   int i, j;
   FILE *f;
   
   //Leitura dos parâmetros do arquivo: tamanho da matriz, elemento do vetor
   //a ser mostrado ao final, erro máximo permitido e número máximo de iterações.
   f = fopen(argv[1], "r");
   fscanf(f, "%d", &n);
   fscanf(f, "%d", &row);
   fscanf(f, "%f", &max_error);
   fscanf(f, "%d", &max_iterations);

   //Dimensões dos blocos e da grid. Os blocos são dispostos linearmente, assim
   //como a grid. Porém, como o tamanho dos blocos é mais limitado que o da grid,
   //cada bloco tem o tamanho máximo BLOCK_SIZE e a grid tem tamanho igual ao
   //lado da matriz, dividindo o processamento de cada linha para um bloco.
   dim3 dimBlock(BLOCK_SIZE, 1, 1);
   dim3 dimGrid(n, 1, 1);
   
   
   //Alocação de memória do host
   hA = (float*) malloc(n * n * sizeof(float*));
   hb = (float*) malloc(n * sizeof(float));
   hx = (float*) malloc(n * sizeof(float));
   herrors = (float*) malloc(n * sizeof(float));
   
   //Alocação de memória da placa de vídeo
   cudaMalloc((void**) &dA, n * n * sizeof(float));
   cudaMalloc((void**) &db, n * sizeof(float));
   cudaMalloc((void**) &dB, n * n * sizeof(float));
   cudaMalloc((void**) &dg, n * sizeof(float));
   cudaMalloc((void**) &dx[0], n * sizeof(float));
   cudaMalloc((void**) &dx[1], n * sizeof(float));
   cudaMalloc((void**) &derrors, n * sizeof(float));
 
   //Leitura dos elementos da matriz e vetor
   for(i = 0; i < n; i++)
   {
      for(j = 0; j < n; j++)
         fscanf(f, "%f", &hA[i * n + j]);
   }
   for(i = 0; i < n; i++)
      fscanf(f, "%f", &hb[i]);

   //Cópia dos elementos lidos do arquivo para a placa de vídeo
   cudaMemcpy(dA, hA, n * n *sizeof(float), cudaMemcpyHostToDevice);
   cudaMemcpy(db, hb, n * sizeof(float), cudaMemcpyHostToDevice);  

   //Cálculo de B e g, inicialização de x
   setup<<<dimGrid, dimBlock>>>(dA, db, dB, dg, dx[0], n);
   
   error = max_error;
   //Cada iteração realiza um update e calcula os erros. Realiza-se então uma
   //redução do vetor de erros em CPU, obtendo-se o maior dentre eles. O loop
   //repete-se enquanto o erro calculado for maior que o permitido e o número
   //de iterações não tenha atingido o máximo.
   for(iterations = 0; iterations < max_iterations && error >= max_error; iterations++)
   {
      update<<<dimGrid, dimBlock>>>(dB, dg, dx[(iterations + 1) % 2], dx[iterations % 2], n);
      errors<<<dimGrid, dimBlock>>>(dA, dx[(iterations + 1) % 2], db, derrors, n);
      cudaMemcpy(herrors, derrors, n * sizeof(float), cudaMemcpyDeviceToHost); 
      error = 0;
      for(i = 0; i < n; i++)
      {
         if(herrors[i] > 0 && herrors[i] > error)
            error = herrors[i];
         if(herrors[i] < 0 && -herrors[i] > error)
            error = -herrors[i];
      }
   }
   /*
   cudaMemcpy(hx, dx[iterations % 2], n * sizeof(float), cudaMemcpyDeviceToHost);
   for(i = 0; i < n; i++)
         printf("x[%d] = %f\n", i, hx[i]);
   */
   printf("Iteracoes: %d\n", iterations);
   printf("Linha %d: calculado = %f, esperado = %f.\n\n", row, hb[row - 1] - herrors[row - 1], hb[row - 1]);
   
   //Libera memória da placa de vídeo
   cudaFree(dA);
   cudaFree(db);
   cudaFree(dB);
   cudaFree(dg);
   cudaFree(dx[0]);
   cudaFree(dx[1]);
   cudaFree(derrors);
   
   //Término da contagem de tempo
   gettimeofday(&tf, NULL);
   sec = tf.tv_sec - ti.tv_sec;
   if(tf.tv_usec < ti.tv_usec)
      sec--;
   usec = (tf.tv_usec - ti.tv_usec + 1000000) % 1000000;
   printf("Tempo decorrido: %ld.%06ld s.\n", sec, usec);
   return 0;
}
